#include "Wavelet.h"
#include <QDebug>
#include <cmath>

Wavelet::Wavelet(const QString& wname) {
    getFilters(wname, Lo_D, Hi_D, Lo_R, Hi_R);
}

void Wavelet::getFilters(const QString& wname, QVector<double>& Lo_D, QVector<double>& Hi_D, QVector<double>& Lo_R, QVector<double>& Hi_R) {
    if (wname == "db1") {
        Lo_D = {0.7071067811865476, 0.7071067811865476};
        Hi_D = {-0.7071067811865476, 0.7071067811865476};
        Lo_R = {0.7071067811865476, 0.7071067811865476};
        Hi_R = {0.7071067811865476, -0.7071067811865476};
    } else {
        qWarning() << "Unsupported wavelet name:" << wname;
    }
}

void Wavelet::dwt(const QVector<double>& x, const QVector<double>& Lo_D, const QVector<double>& Hi_D, QVector<double>& approx, QVector<double>& detail) {
    int len = x.size();
    int filterLen = Lo_D.size();
    int approxLen = (len + filterLen - 1) / 2;
    int detailLen = (len + filterLen - 1) / 2;

    approx.resize(approxLen);
    detail.resize(detailLen);

    // 近似系数（低通滤波）
    for (int i = 0; i < approxLen; i++) {
        double sum = 0.0;
        for (int j = 0; j < filterLen; j++) {
            int idx = 2 * i - j;
            if (idx >= 0 && idx < len) {
                sum += Lo_D[j] * x[idx];
            }
        }
        approx[i] = sum;
    }

    // 细节系数（高通滤波）
    for (int i = 0; i < detailLen; i++) {
        double sum = 0.0;
        for (int j = 0; j < filterLen; j++) {
            int idx = 2 * i - j;
            if (idx >= 0 && idx < len) {
                sum += Hi_D[j] * x[idx];
            }
        }
        detail[i] = sum;
    }
}

void Wavelet::idwt(const QVector<double>& approx, const QVector<double>& detail, const QVector<double>& Lo_R, const QVector<double>& Hi_R, QVector<double>& x) {
    int len = approx.size()  + detail.size();
    int filterLen = Lo_R.size();
    x.resize(len);

    // 重构信号
    for (int i = 0; i < len; i++) {
        double sum = 0.0;
        for (int j = 0; j < filterLen; j++) {
            int idx = (i - j) / 2;
            if (idx >= 0 && idx < approx.size()  && (i - j) % 2 == 0) {
                sum += Lo_R[j] * approx[idx];
            }
            if (idx >= 0 && idx < detail.size()  && (i - j) % 2 == 0) {
                sum += Hi_R[j] * detail[idx];
            }
        }
        x[i] = sum;
    }
}

void Wavelet::wavedec(const QVector<double>& x, int n, QVector<double>& c, QVector<int>& l) {
    c.clear();
    l.resize(n  + 2);
    l[n + 1] = x.size();

    QVector<double> approx = x;
    QVector<double> detail;

    for (int k = 1; k <= n; k++) {
        dwt(approx, Lo_D, Hi_D, approx, detail);
        c = detail + c;
        l[n + 1 - k] = detail.size();
    }

    c = approx + c;
    l[0] = approx.size();
}

void Wavelet::waverec(const QVector<double>& c, const QVector<int>& l, QVector<double>& x) {
    QVector<double> approx = c.mid(0,  l[0]);
    QVector<double> detail;

    for (int k = 1; k < l.size()  - 1; k++) {
        detail = c.mid(l[0],  l[k]);
        idwt(approx, detail, Lo_R, Hi_R, approx);
    }

    x = approx;
}
